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ABSTRACT 

We study the population statistics of the surviving subhaloes of ACDM dark matter 
haloes using a set of very high resolution TV-body simulations. These include both 
simulations of representative regions of the Universe and ultra-high resolution resim- 
ulations of individual dark matter haloes. We find that more massive haloes tend to 
have a larger mass fraction in subhaloes. For example, cluster size haloes typically 
have 7.5 percent of the mass within -R200 m substructures of fractional mass larger 
than 10 -5 , which is 25 percent higher than galactic haloes. There is, however, a large 
variance in the subhalo mass fraction from halo to halo, whereas the subhalo abun- 
dance shows much higher regularity. For dark matter haloes of fixed mass, the subhalo 
abundance decreases by 30 percent between redshift 2 and 0. The subhalo abundance 
function correlates with the host halo concentration parameter and formation redshift. 
However, the intrinsic scatter is not significantly reduced for narrow ranges of con- 
centration parameter or formation redshift, showing that they are not the dominant 
parameters that determine the subhalo abundance in a halo. 
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1 INTRODUCTION 

In the standard ACDM model of cosmogony, cold dark mat- 
ter haloes form from primordial Gaussian fluctuations by 
accretion of diffuse dark matter and mergers with other 
haloes. During this hierarchical process, accreted haloes of- 
ten survive as self-bound substructures or subhaloes or- 
biting around the main system (e. g 



2004 



2005 



Tormen et al. 199 



Ghigna et all l2000l: ISpringel et al.l l200ll: |Pe Lucia et al 



,Gao et al.l l2004al; IDiemand et al.l]2004l; iReed et~, 
iDiemand et al.ll2007l ; ISpringel et aljl2008af ). The most 
massive subhaloes are likely hosts of luminous satellite 
galaxies. 

The evolution of subhaloes is determined by the inter- 
play of phenomena such as accretion, mergers, dynamical 
friction and tidal stripping. Thus, subhalo properties are 
most directly and accurately investigated using simulations. 
So long as one is exclusively interested in evolution of the 
dark matter, this is a purely gravitational problem that is 
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ideally suited to N-body techniques. Indeed, most work on 
subhaloes has exploited this approach. However, since sub- 
haloes represent a small fraction of the total mass of the 
halo, most studies have been based, for computational re- 



of individual haloes (e.g.jGhigna et al.ll2000: Springel et al. 


2001 


; Stoehr et al.l 20031 


Gao et al. 2004al: IDiemand et al. 


2004 


. 120071: Springel et al. 


2008a|). 



Statistical studies are much more demanding because 
of the need to simulate large volumes with a large dy- 
namic range and, to date, have focused on subhaloes in rel- 

ative massive haloes (e.g 
■ 1 1 -1 1—"- 

2004 
I2OOI 



De Lucia et al. 20041: iGao et al 



Kravtsov et al] 120041 : IShaw et alj l2006l ; lElahi et al 
lAngulo et all 120091 ). Inspired by, and calibrated 
from, the results of simulations, semi-analytic models have 
been developed to c alculate certa i n properties of sub - 
halo populations (e.g., [B enson 20 051: iTavlor fc Babu l 200 



den Bosch et al.ll2005l : IZentner et al.ll2005l : lGiocoli et al 
2008al lbl: iLi fc Mo 1120091 ). 

In this short paper, we exploit the unprecedented statis- 
tical power and high resolution of the Millennium-II Simula- 
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tion (|Bovl an-Kolchin et al ]|2009l ), together with other recent 
cosmological simulations and a suite of very high resolution 
resimu lations of individual galaxy haloes (jSpringel et al.l 
2008a) to study the abundance of subhaloes with much bet- 
ter statistics and over a significantly larger dynamic range 
than has been possible in the past. We address three spe- 
cific questions concerning the subhalo mass function: (i) how 
does it depend on the mass of the parent halo? (ii) how does 
it depend on redshift, for a fixed halo mass? (iii) how does it 
depend, at fixed mass and redshift, on the properties of the 
halo, specifically the concentration and formation redshift? 
Some of these questions have been addressed in some of the 
earlier studies mentioned above. Our results extend these 
studies into a regime of subhalo mass not previously probed 
and with much greater statistical power. 



Table 1. Simulation parameters: (1) name of the simulation; (2) 
the side of the simulated cubic volume; (3) the number of particles 
in the simulation; (4) the mass per particle; (5) the gravitational 
softening length. The initial conditions for hMS and MS-II simu- 
lations have matching phases. 



Name 


Box size 
[h- 1 Mpc] 


N p 


r»M ] 


e 

[h'^-kpc] 


MS 


500 


2160 3 


8.6 x 10 s 


5.0 


hMS 


100 


900 3 


9.5 x 10 7 


2.4 


MS-II 


100 


2160 3 


6.9 x 10 6 


1.0 


Aquarius 
(level 2) 






~ 1 x 10 4 


0.048 



2 NUMERICAL SIMULATIONS 

The main advance in this paper comes from analysing 
two sets of very recent, high resolution simulations. One 
is the 10-billion particle M illennium-II Simulation (MS-II, 
iBovlan-Kolchin et al.ll2009l ) of a 100/i _1 Mpc cubic volume. 
The oth er one is the set of 6 r esimulations in the Aquarius 
Project |Springel et al. I l2008al lbh. each of which follow the 
evolution of a galaxy size halo with more than 10 8 parti- 
cles inside the virial radius. All these simulations represent 
a significant improvement over the previous generation of 
N-body simulations. We supplement these data with older 
simulations, t he original 10-billion particle Millennium Sim- 
ulation (MS), ISpringel et al.l I2005T) of a 500fe~* Mpc cubic 
volume and the hMS simulation (|Gao et al.ll2008l ). which has 
matching phases in the initial conditions to those of the MS- 
II, but with about 14 times lower mass resolution and 2.4 
times lower spatial resolution. 

All the simulations assume the same values of the cos- 
mological parameters: mean matter density, Q m = 0.25; cos- 
mological constant term, Qa = 0.75; Hubble parameter (in 
units of 100 kms -1 Mpc -1 ), h = 0.73; fluctuation ampli- 
tude, ag — 0.9; and primordial spectrum power-law index, 
n — 1 ; These values are c onsistent with the first year WMAP 
data ijSpergel et a\fo o03) but differ at about the 2 — a level 
from more recent determinations. This small off-set is of 
no consequence for the topics addressed in this p aper. All 
the s imulations were run with either the Gadget-2 (|Springell 
120051 ) or the more recent, Gadget-3, code used for the Aquar- 
ius Project. 

The parameters of the simulations are summarized in 
Table [T] Note that the MS and the MS-II differ in mass res- 
olution by a factor of 125, while the MS and the hMS differ 
by a factor of 9, allowing numerical convergence tests to be 
carried out. 

Dark matter hal oes were found u sing the friends- 
of- friends algorithm l|Davis et al.l 1 19851 ). and their self- 
bound sublialoes_werejdentified using the SUBFIND algo- 
rithm jSpringel et al .l200ll ). Haloes and their subhaloes were 
tracked across the simulation outputs and linked together in 
merger trees. In what follows, we restrict our attention to 
subhaloes contained within the virial radius, -R200, of their 
parent host halo, defined as the radius within which the 
enclosed mass density is 200 t i mes t he critical value. We 
note that IBovlan-Kolchin et al.l l|2010l ) used a different def- 
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Figure 1. The mean cumulative mass function of subhaloes in 
hosts with masses in the range [1,3] X 10 14 /i~ 1 Mq, as a func- 
tion of subhalo mass fraction, M su ^/M200- The y-axis has been 
multiplied by the mass fraction in order to expand the dynamic 
range. The 3 solid lines show the averaged subhalo mass function 
for samples of 12 haloes in the MS-II (red) and the hMS (blue) 
and of 1676 haloes in the MS (black), as indicated in the legend. 
An analytic fit to the MS data using eqn. JTJ is shown as a dotted 
line. The error bars on selected points show the rms scatter about 
the mean for the MS and MS-II haloes. 

inition of the virial radius, base d on the spherical collapse 
model l|Eke. Cole fc Frenklll996l ). 



3 RESULTS 

We start by considering the cumulative subhalo mass func- 
tions for samples of haloes with mass in the range [1,3] x 
1O 14 /i _1 M in the MS-II, hMS and MS simulations at z = 0. 
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The first two samples comprise the same 12 haloes at differ- 
ent resolutions, while the third contains 1676 haloes. 

The mean cumulative mass functions of subhaloes for 
all three samples are plotted in Fig. [JJ Here, and below, we 
show cumulative mass functions multiplied by M sub /M2oo 
in order to take out the dominant mass dependence and 
make the differences be tween curves more apparent. Follow- 
ing previous work (e .g., Stoehr et all2003l; 



iDiemand et al.1 [20071 ; ISpringel et al. 12008a ). we restrict the 
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samples to include only subhaloes with 100 particles or more 
in order to avoid numerical effects. There is good conver- 
gence between the results from the MS-II and the hMS and, 
within the errors, there is good agreement with the cumula- 
tive mass function of the MS. The error bars show the rms 
scatter about the mean for the MS and MS-II haloes. (They 
are plotted at a selection of well-separated points to reduce 
the effect of correlations between adjacent bins.) This scat- 
ter, which reflects the varied formation histories of haloes, 
is substantial, particul arly for large values of M su b /M2oo. 

In a recent study, iBovlan-Kolchin et al.l l|2010T ) demon- 
strated that the subhalo occupation distribution is well de- 
scribed by a negative binomial, independently of host halo 
mass. In this case, the scatter at the high mass end of the 
subhalo mass function is well described by a Poisson dis- 
tribution, but at the low mass end, m aub /M2oo < 10 -4 , in- 
trinsic scatter of fractional amplitude a ~ 0.15 dominates. 
The authors further found that the cumulative subhalo mass 
function of galactic size haloes can be well fitted with the 
following form: 



N( >f i 



/M 2 



exp 



^CUt 



(1) 



The dashed line in Fig. [TJ is a fit to the mean subhalo 
cumulative mass function for the cluster size haloes from the 
MS using eqn. ([TJ with parameters Jii — 0.01, fi cu t = 0.1, 
b — 1.2 and power-law index a = —0.94. Clearly, the fit 
is very good. The best-fit parameters differ sli ghtly from 
the values given by IBovlan-Kolchin et al.l l|2010h because a 
different definition of virial radius has been used here and 
the range of halo masses is different. 



3.1 Dependence on host halo mass 

iGao et ahl (|2004aT ) noted that the abundance of relatively 
massive subhaloes (M su b/M2oo > 0.001) increases system- 
atically with host halo mass. This trend reflects the fact 
that in the CDM cosmology more massive haloes are, in the 
mean, both less centrally concentrated and younger than 
less massive haloes. Thus, they exert weaker tidal forces and 
have had less time to disrupt their substructure. This re- 
sult has been confirmed by subsequent numerical and semi- 
analytical studies JZentner et al .1120051: Ivan den Bosch et al.l 
l2005l : IShaw et al.ll2006l : lGiocoli et al.ll2008bh . 

The much larger samples of well resolved haloes in our 
set of simulations allows us to explore the trend with host 
halo mass to much smaller mass ratios than have been 
resolved in earlier simulations which, f or the most part , 
have focused on cl uster size haloes (e.g. IGao et al 2004a; 
IZentner et al.|[2005l ). Combining data from the Aquarius and 
MS-II simulations, we have a sample of haloes spanning the 
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Figure 2. The dependence of the subhalo mass function on 
the host halo mass. The solid lines show the averaged cumu- 
lative subhalo mass functions for three intervals of host halo 
mass: [1,3] X 10 14 r'Mg (black), [1,3] X IO^/i^Mq (red) and 
[1,3] X 10 12 h~ 1 Mq (blue). Haloes in the first mass range come 
from the MS and MS-II, while those in the latter two ranges come 
from the MS-II. The error bars on selected points show the er- 
ror on the mean for the three mass ranges indicated. The filled 
squares show the mean of the cumulative subhalo mass functions 
of the 6 Aquarius haloes. The dotted lines show fits to eqn. ((TJ; 
the fit parameters are listed in eqn. f2l . 



range ~ 10 12 — 3 x 10 14 /i _1 Mq, each resolved with at least 
10 5 particles within the virial radius, and often with 10 to 
1000 times more. This represents an improvemen t by a fac- 
tor of 1 00-1000 over the simulations analyzed bv lGao et al.l 
(|2004aT ). 

Fig. [2] shows cumulative subhalo mass functions for 
host haloes of different mass. The magenta squares and 
the solid lines respectively represent galactic size haloes 
([1,3] x 10 12 /i _1 M Q ) from Aquarius and MS-II, the blue 
lines represent group size haloes ([1,3] x 10 13 /i _1 Mq) from 
the MS-II and the black lines represent cluster size haloes 
([1,3] x 10 14 /i _1 M©) from the MS-II and the MS. The cor- 
responding best fits of the form of eqn. ([TJ are overplotted 
as dotted lines. Clearly, eqn. {TJ provides an excellent fit to 
the cumulative subhalo mass function for all three ranges of 
host halo mass, albeit with varying best-fit parameters. We 
fixed the subhalo function slope index to a — —0.94 when 
fitting all three curves. Of the remaining three free parame- 
ters in eqn. {TJ, ~\±l determines the overall amplitude of the 
subhalo abundance function and is largely independent of 
the parameters /i cu t and b which combine to determine the 
shape of the cut-off at the high mass end. We find that nei- 
ther fi cut nor b can be well constrained on their own, but 
are largely degenerate. Choosing fj, cut in the range [0.04, 0.1] 
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Figure 3. The mass fraction in substructures in haloes of differ- 
ent mass as a function of the normalised subhalo mass. The lines 
show the medians of the cumulative subhalo mass fractions for 
three ranges of host halo mass in the MS-II simulation: [1, 3] X 10 14 
/i -1 M Q (black solid lines), [1,3] X 10 13 /i _1 M Q (red solid lines) 
and [1,3] X 10 12 /i —1 MQ (blue solid lines). The magenta squares 
show the averaged subhalo mass fraction in the 6 Aquarius haloes. 
The dashed lines show the 20 and 80% of the corresponding dis- 
tribution. The error bars on selected points show the error on the 
median. 



Figure 4. The dependence of the subhalo mass function on red- 
shift. The lines show the averaged cumulative subhalo mass func- 
tions for MS-II haloes in the mass range [1 — 3] X 10 13 /i — 'Mq 
at redshift 2 = (black), z = 0.5 (red), 2 = 1 (blue) and 2 = 2 
(green). The filled squares show results for the Aquarius haloes at 
2 = 1 (large red) and 2 = (magenta) . The error bars on selected 
points show the error on the mean. 



or b in the range [0.8, 1.5] gives fits which match the overall 
subhalo halo mass function reasonably well. To break the 
degeneracy between ^t cut and b in our fits we set b = 1.2 
and fitted eqn. fl]) by varying just /ii and /i cu t. The best-fit 
values for /Ii and ^ cu t are: 

M 200 € [1, 3] x 10 14 h" 1 M : jE = 0.0110, /w = 0.10. 
M200 G [1, 3] x 10 13 h" 1 M Q : /IT = 0.0092, Mcut = 0.07. 
M200 G [1, 3] x 10 12 h' 1 M Q : /IT = 0.0085, /x cut = 0.08. (2) 

By examining the subhalo mass functions of individual 
halos, we have found that the main scatter is in the nor- 
malisation rather than the shape; halos with a higher than 
average abundance of low-mass subhalos also tend to have a 
higher than average abundance of high-mass halos and vice 
versa. 

At a given M su b/M2oo, there is a weak trend in the 
abundance of subhaloes with host halo mass. In the region 
of overlap in M su b/M2oo, the cluster haloes have a mass frac- 
tion of substructures that is about 25% higher than in the 
galaxy haloes. The mean of the group sample is intermedi- 
ate between these two. The 1 5% differen c e betw een group 
and cluster halo identified by iGao et a l. (2004a) is visible 
at the smaller mass ratios plotted, M su b/M2oo ~ 10 -4 . The 
mass dependence of the subhalo abundance with host halo 



mass in our simulations is much weaker than a theoretical 
expectation based on semianalytical modelling of the sub- 
halo population (e.g.. Ivan den Bosch et al.ll2005h . 

The mass fractions in subhaloes as a function of relative 
mass for host haloes of different size are shown in Fig. [3] The 
curves are very steep at high values reflecting the fact that 
most of a typical halo's mass in subhaloes is contributed 
by a relati vely small number of v ery massive subhaloes. As 
shown by ISpringel et al.l l|2008al ) for the Aquarius haloes, 
these massive subhaloes tend to be in the outer parts of their 
parent halo. The subhalo mass fraction above a given value 
of M su b/Af2uo depends on the mass of the parent halo. Typ- 
ically, cluster haloes contain about 25% more subhalo mass 
than galaxy haloes. The Aquarius galactic haloes resolve 
substructure with relative mass as small as 10 -6 . Subhaloes 
above this mass contain about 7% of the total halo mass 
within i?200- For clusters, the MS-II resolves subhaloes with 
relative mass just below 10 -5 . Subhaloes more massive than 
this already amount to about 7% of the mass within i?200- 
Note, however, that there is a rather large variance from halo 
to halo in the subhalo mass fraction above a given value of 
M su b/M2oo, particularly at the high mass end, where much 
of the subhalo mass is typically contributed by one or two 
objects. 
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3.2 Redshift dependence 

We expect haloes of a given mass to contain more subhaloes 
at earlier times because the earlier counterparts are both 
less concentrated and dy namically younger. This trend was 
seen in the simulations of iGao et al,l ((2004a) . To investigate 
the redshift dependence of the subhalo mass function in our 
simulations, we use a sample of group haloes of mass [1,3] x 
10 13 /i _1 M Q in the MS-II. We also use the Aquarius galaxy 
haloes at z = and their main progenitors at z = 1, the 
mean mass of these halo is M200 = 1.19 x W 12 h~ 1 MQ at 
z = and M 20 o = 6.73 x lO n /i _1 M at z = 1. 

The results are displayed in Fig [4] In the MS-II sam- 
ple, there are a total of 219, 204, 169 and 41 haloes at the 
four redshifts shown, 2 = 0, 0.5, 1, 2. Clearly, the subhalo 
mass function evolves with redshift weakly, but systemat- 
ically, in the expected sense: at high redshift haloes of a 
given mass contain more subhaloes. For example, the sub- 
halo abundance in groups at z = is typically 18% lower 
than at z — 0.5, 25% lower than at z = 1, and 30% lower 
than at z — 2. The subhalo mass function of galaxy haloes 
evolves slightly more slowly than this, reflecting the earlier 
formation epoch of galaxy haloes. In this case the abundance 
is only 15% higher at z = 1 than at z = 0. We note that 
using eqn. (QJ with a and b fixed to —0.94 and 1.2, respec- 
tively, fits the data shown at each redshift very well. The 
corresponding parameters are: 

z = 0.0: JTi = 0.0092, fi cut = 0.07 
z = 0.5: JTi = 0.0118, /Ucut = 0.06 
z = 1.0: Jii = 0.0130, A*cut = 0.05 

z = 2.0: JTi = 0.0140, ^cut = 0.02 (3) 

For the galaxy haloes the trend persists down to the smallest 
subhalo masses resolved in the set of six Aquarius simula- 
tions, M sub /M 20 o = 10" 6 . 

3.3 Dependence on host halo properties 

We now examine how the subhalo mass function depends 
on two basic properties of the parent halo: the concentra- 
tion parameter and the formation redshift. We also consider 
whether any such dependence contributes to the scatter in 
the subha lo mass functi o n see n in previous figures, as sug- 
gested bv lZentner et all (|2005l ). Earlier studies have shown 
that, at least for relatively massive haloes, the subhalo abun- 
dance does depend on the concentration and formation red- 
shift of the parent ha lo (|Gao et al.l l2004al : IZentner et al.l 
l2005l : IShaw et al.ll2006T ). With the MS-II, we can re-examine 
the dependence of subhalo abundance on the properties of 
the parent halo with much better statistics than was possi- 
ble before and also explore the low mass end of the subhalo 
mass distribution. 

We select 219 haloes in the MS-II in the mass range 
[1,3] x 1O 13 /i _1 M0, and subdivide this sample into three 
equal size subsamples ranked according to concentration 
parameter or formation redshift. We evaluate the concen- 
tr ation parameter o f a halo as Vmax/V^oo (as was done 
111 IGao et al l [200481 ) and we define its formation redshift 
as the time when half the mass was asse mbled. The result s 
shown in Fig. [5] confirm the conclusion of lGao et al.l (|2004al ) 
that the subhalo abundance decreases with increasing par- 



ent halo concentration and formation redshift. At the low 
subhalo mass end, the third least concentrated and most 
recently formed haloes have about 25% more substructures 
than the most concentrated and earliest forming third. At 
the high subhalo mass end, Msub/Afeio > 0.01, the differ- 
ence between the two extreme thirds is about a factor of 2, 
significantly larger than at the lower subhalo mass end. 

The variance around the mean of the subhalo mass func- 
tion for each subsample and for the sample as a whole are 
listed in Table [21 in units of the mean value. The size of the 
trends seen in Fig. [S]is comparable to the scatter in the re- 
lations. Nevertheless, the dependence of subhalo abundance 
on the concentration and formation time of the parent halo, 
at a fixed host halo mass, is much stronger than the depen- 
dence on host halo mass seen in Fig. [2] 

The data presented in Table [2] show that the scatter in 
subhalo abundance among haloes chosen to have a narrow 
range of concentration parameters or formation times is not 
significantly smaller than the scatter for the halo population 
as a whole. Thus, neither concentration nor formation time 
appear to be the primary halo property responsible for the 
observed scatter in the subhalo mass function. 



4 CONCLUSION 

We have used a new set of very large cosmological TV-body 
simulations to investigate the statistics of subhalo abun- 
dance in ACDM dark matter haloes. Our results may be 
summarized as follows: 

(i) The subhalo abundance function of dark matter haloes 
can be well fitted w it h the functional form proposed by 
iBovlan-Kolchin et al.l (|2010h . which is a power-law at the 
low subhalo mass end and has an exponential cut-off at the 
high mass end, independently of host halo mass and redshift. 

(ii) The subhalo abundance function depends weakly on 
host halo mass. The difference between typical cluster and 
galaxy size haloes is about 25 percent at the lower subhalo 
mass end, substantially weaker than expected from previous 
studies. The scatter in the subhalo abundance function of 
different haloes is larger than the systematic trend. 

(iii) For a given mass halo, the subhalo abundance evolves 
systematically with redshift. The evolution is largest at 
lower redshift (~ 18 percent between z — and z = 0.5) 
and becomes very small at high redshift (a few percent be- 
tween z = 1 and z — 2). Over the range z — — 2, the 
substructure abundance increases by 30% for haloes corre- 
sponding to bright galaxies and poor galaxy groups. 

(iv) At fixed mass, dark matter haloes with higher con- 
centration parameters and earlier formation redshifts con- 
tain fewer subhaloes. The difference is a factor of 25 percent 
between the top and bottom thirds of the population ranked 
by concentration parameter or formation redshift. However, 
this dependence on concentration and formation redshift is 
not enough to explain the scatter in the subhalo mass func- 
tion, suggesting that these are not the dominant parameters 
determining the subhalo mass fraction. 
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Figure 5. The dependence of the subhalo mass function on the properties of the parent halo. The left panel shows the dependence 
on the concentration parameter of the halo and the right panel on the formation redshift. The lines show the mean of the cumulative 
subhalo mass functions for 219 MS-II haloes in the mass range [1 — 3] X 1O 13 /i _1 M0. In the two panels, the red solid lines show the 
averaged subhalo abundance function of the top third of the sample, that with the highest concentration parameter (left), and earliest 
formation redshift (right). The blue lines correspond to the intermediate third of the sample and the black lines to the bottom third, 
with the lowest concentration parameter and the latest formation redshift. The error bars on the selected points show the error on the 
mean. 



Table 2. The standard deviation from the mean of the subhalo mass function at a given fractional mass, m su b / M200 = 0.01, 0.001, 0.0001, 
in units of the mean. 



fi high c intermediate c low c high z intermediate z low z whole sample 



0.01 1.52 1.23 0.67 1.5 1.01 0.81 1.08 

0.001 0.47 0.44 0.32 0.48 0.43 0.33 0.48 

0.0001 0.21 0.21 0.18 0.20 0.20 0.20 0.22 
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